The data presented here was sampled from movebank.org and was originally published here:
Carrasco-Harris MF, Bowman D, Reichling S, Cole JA. 2020. Spatial ecology of copperhead snakes (Agkistrodon contortrix) in response to urban park trails. J Urban Ecol. 6(1):juaa007. https://doi.org/10.1093/jue/juaa007
In their study, the authors measured the spatial activity of copperheads (Agkistrodon contortrix) around pedestrian trails in an urban forest (Overton Park, Memphis, TN) using radio telemetry. They found that sex and season, not distance to the nearest trail, affected the distance that snakes moved and concluded that copperheads may be tolerant of low-level human disturbances.
Import Dataset with
read.csv
snake.data <- read.csv("data/snakes.csv")
snake <- snake.data %>%
filter(individual.local.identifier == c("ID_4241","ID_4270","ID_7270","ID_4266"))longer object length is not a multiple of shorter object length
I decided to subset for just 4 individuals that I thought gave a good representation of the data. These snakes each primarily resided in different sections of the park and displayed variation in the sizes of their respective home ranges.
The points appear to be have a similar spread for each individual without any clear outliers among them (at least to my eyes).
qaqc_plot <- ggplot() + geom_point(data=snake,
aes(utm.easting,utm.northing,
color=individual.local.identifier)) +
labs(x="Easting", y="Northing") +
guides(color=guide_legend("Identifier"))
ggplotly(qaqc_plot)lapply(split(snake, snake$individual.local.identifier),
function(x)write.csv(x, file = paste(x$individual.local.identifier[1],".csv", sep = ""), row.names = FALSE))$ID_4241
NULL
$ID_4266
NULL
$ID_4270
NULL
$ID_7270
NULL
utm_points <- cbind(snake$utm.easting, snake$utm.northing)
utm_locations <- SpatialPoints(utm_points,
proj4string=CRS("+proj=utm +zone=16 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster2 <- openmap(c(max(proj_lat.lon$y)+0.003, min(proj_lat.lon$x)-0.003),
c(min(proj_lat.lon$y)-0.003, max(proj_lat.lon$x)+0.003),
type = "bing")
raster <- openmap(c(max(proj_lat.lon$y)+0.01, min(proj_lat.lon$x)-0.01),
c(min(proj_lat.lon$y)-0.01, max(proj_lat.lon$x)+0.01),
type = "bing")
raster_utm <- openproj(raster,
projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")
raster2_utm <- openproj(raster2,
projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")I decided to make a second raster that still takes into account the location of the points when choosing the maximum and minimum coordinates of the resulting map, but with a smaller window to get us slightly more zoom over our park (which represents such a small area).
Here’s the map made with the raster inputs from Dr. Gentry’s code:
autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
theme(legend.position="right") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
geom_point(data=snake, aes(utm.easting,utm.northing,
color=individual.local.identifier), size = 1, alpha = 0.8) +
theme(axis.title = element_text(face="bold")) + labs(x="Easting",
y="Northing") + guides(color=guide_legend("Identifier"))And with a constrained window:
autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() +
theme(legend.position="right") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
geom_point(data=snake, aes(utm.easting,utm.northing,
color=individual.local.identifier), size = 1, alpha = 0.8) +
theme(axis.title = element_text(face="bold")) + labs(x="Easting",
y="Northing") + guides(color=guide_legend("Identifier"))For some reason I was having an issue with the permit status (?) using the pbapply package. For whatever reason and completely mysteriously to me, the first three lines of code in this chunk fixed the issue for me
library(maptools)
if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit() #This resolved the following error message from pbapply: "Error: isTRUE(gpclibPermitStatus()) is not TRUE"[1] TRUE
mcp_raster <- function(filename){
data <- read.csv(file = filename)
x <- as.data.frame(data$utm.easting)
y <- as.data.frame(data$utm.northing)
xy <- c(x,y)
data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
xy <- SpatialPoints(data.proj@coords)
mcp.out <- mcp(xy, percent=100, unout="ha")
mcp.points <- cbind((data.frame(xy)),data$individual.local.identifier)
colnames(mcp.points) <- c("x","y", "identifier")
mcp.poly <- fortify(mcp.out, region = "id")
units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85, y=0.95,
gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
mcp.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), alpha=0.5, fill = "red") +
geom_point(data=mcp.points, aes(x=x, y=y)) +
labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
annotation_custom(units)
mcp.plot
}
pblapply(files, mcp_raster)
| | 0 % ~calculating
|================ | 25% ~01s
|=============================== | 50% ~00s
|============================================== | 75% ~00s
|=============================================================| 100% elapsed=01s
[[1]]
[[2]]
[[3]]
[[4]]
kde_raster <- function(filename){
data <- read.csv(file = filename)
x <- as.data.frame(data$utm.easting)
y <- as.data.frame(data$utm.northing)
xy <- c(x,y)
data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
xy <- SpatialPoints(data.proj@coords)
kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
ver <- getverticeshr(kde, 95)
kde.points <- cbind((data.frame(data.proj@coords)),data$individual.local.identifier)
colnames(kde.points) <- c("x","y","identifier")
kde.poly <- fortify(ver, region = "id")
units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85, y=0.95,
gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
kde.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.5, fill = "red") +
geom_point(data=kde.points, aes(x=x, y=y)) +
labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
annotation_custom(units)
kde.plot
}
pblapply(files, kde_raster)
| | 0 % ~calculating
|================ | 25% ~02s
|=============================== | 50% ~01s
|============================================== | 75% ~00s
|=============================================================| 100% elapsed=01s
[[1]]
[[2]]
[[3]]
[[4]]
snake_4241 <- read.csv("ID_4241.csv")
date <- as.POSIXct(strptime(as.character(snake_4241$timestamp),"%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"))
snake_4241$date <- date
snake_4241.reloc <- cbind.data.frame(snake_4241$utm.easting, snake_4241$utm.northing,
as.vector(snake_4241$individual.local.identifier),
as.POSIXct(date))
colnames(snake_4241.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(snake_4241.reloc, date=date, id="snake_4241")
sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
snake_4241.traj <- kernelbb(trajectory, sig1 = .7908, sig2 = 58, grid = 100)
bb_ver <- getverticeshr(snake_4241.traj, 58)
bb_poly <- fortify(bb_ver, region = "id",
proj4string = CRS("+proj=utm +zone=16+
datum=WGS84 +units=m +no_defs"))
colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")
bb_image <- crop(snake_4241.traj, bb_ver,
proj4string = CRS("+proj=utm +zone=16 +
datum=WGS84 +units=m +no_defs"))
bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85, y=0.95,
gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
bb.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
geom_tile(data=bb_image,
aes(x=bb_image@coords[,1], y=bb_image@coords[,2],
fill = bb_image@data$ud)) +
geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "black", fill = NA) +
scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
labs(x="Easting (m)", y="Northing (m)", title="4241") +
theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
bb.plotsnake_4270 <- read.csv("ID_4270.csv")
date <- as.POSIXct(strptime(as.character(snake_4270$timestamp),"%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"))
snake_4270$date <- date
snake_4270.reloc <- cbind.data.frame(snake_4270$utm.easting, snake_4270$utm.northing,
as.vector(snake_4270$individual.local.identifier),
as.POSIXct(date))
colnames(snake_4270.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(snake_4270.reloc, date=date, id="snake_4270")
sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
snake_4270.traj <- kernelbb(trajectory, sig1 = .7908, sig2 = 58, grid = 100)
bb_ver <- getverticeshr(snake_4270.traj, 24.5)
bb_poly <- fortify(bb_ver, region = "id",
proj4string = CRS("+proj=utm +zone=16+
datum=WGS84 +units=m +no_defs"))
colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")
bb_image <- crop(snake_4270.traj, bb_ver,
proj4string = CRS("+proj=utm +zone=16 +
datum=WGS84 +units=m +no_defs"))
bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85, y=0.95,
gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
bb.plot <- autoplot.OpenStreetMap(raster2_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
geom_tile(data=bb_image,
aes(x=bb_image@coords[,1], y=bb_image@coords[,2],
fill = bb_image@data$ud)) +
geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "black", fill = NA) +
scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
labs(x="Easting (m)", y="Northing (m)", title="4270") +
theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
bb.plotInterestingly, even though 4270 had a home range that was fairly comparable to 4241 in size (at least by mcp), adding the timestamp data and generating this ‘heat map’ of activity shows that 4270’s activity was generally concentrated in a much denser area.